Example SEND analysis using PyXEM on DLS Jupyterhub¶
author: Mohsen Danaie
This notebook is intended to be run on Diamond's Jupyter hub. See instructions to connect here
Reference: Ophus, C. (2019). Four-Dimensional Scanning Transmission Electron Microscopy (4D-STEM): From Scanning Nanodiffraction to Ptychography and Beyond. Microscopy and Microanalysis, 25(3), 563-582. doi_link
The 4D-STEM technique refers to a wide range experimental set-up in the transmission electron microscope where the electron beam is rastered across a thin sample and an electron diffraction pattern is recorded post-specimen corresponding to each probe position. This notebook addresses a subgroup of this technique, the so-called Scanning Electron Nanodiffraction (SEND). Here, an electron probe is formed with a small convergence semi-angle (~1 mrad) enabling the diffraction discs from crystalline phases within the sample to be well-resolved, (compare panel b above to a where the probe convergence is large resulting in overlap of diffraction features). SEND can be used for phase, orientation, and strain mapping.
The demonstration in this notebook covers how we can directly load a raw dataset (17 GB) and visualise it using pyxem. We then go through the calibration routine performed on the data from standard sample (gold cross-grating) to be then applied to other experimental data.
Importing the libraries needed¶
%matplotlib notebook
import numpy as np
import hyperspy.api as hs
import pyxem as pxm
import matplotlib.pyplot as plt
import os
import sys
from pyxem.libraries.calibration_library import CalibrationDataLibrary
from pyxem.generators.calibration_generator import CalibrationGenerator
import numpy as np
import matplotlib.pyplot as plt
from skimage import data
from skimage.registration import phase_cross_correlation
from skimage.registration._phase_cross_correlation import _upsampled_dft
from scipy.ndimage import fourier_shift
import math
import warnings
warnings.simplefilter(action='ignore')
plt.rcParams.update({'figure.max_open_warning': 0})
warnings.filterwarnings("ignore")
current_path = os.getcwd()
if os.path.exists(os.path.join(current_path, 'epsic_tools')):
sys.path.append(os.path.join(current_path, 'epsic_tools'))
else:
!git clone https://github.com/ePSIC-DLS/epsic_tools
sys.path.append(os.path.join(current_path, 'epsic_tools'))
import epsic_tools.api as epsic
(a) Example PyXem functionalities¶
Loading raw data¶
Using the following method we can load a raw dataset into the notebook lazily:
data = pxm.load_mib('/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/expands_notebook/raw_mib/Au_Xgrating_200kX_CL20cm_10um_256_1ms.mib')
print(data)
Data visualisation by forming virtual dark-field images¶
For visualisation we need to load into memory:
data.compute()
data.plot(vmin=0.,vmax=10)
The green circular region of interest above (this appears after running the cell below), defines a virtual aperture. The sum intensity of this region then is displayed as a virtual dark-field image below. By changing the size and location of the aperture size, the VDF image is updated interactively.
roi = hs.roi.CircleROI(cx=273, cy=319, r=10, r_inner=0)
data.plot_integrated_intensity(roi)
(b) Calibrations¶
Loading files¶
Loading files from the starndard cross-grating gold sample for the calibrations. Note that these files were converted from raw to hdf5 / hspy and tiff formats.
base_folder = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/expands_notebook/'
cals_path = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/expands_notebook/au_xgrating_cals'
mask_path = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/expands_notebook/mask.hspy'
hdf_files = []
hspy_files = []
tiff_files = []
for dirname, dirnames, filenames in os.walk(cals_path):
for filename in filenames:
if os.path.splitext(filename)[1] == '.hdf5':
hdf_files.append(os.path.join(dirname, filename))
if os.path.splitext(filename)[1] == '.hspy':
hspy_files.append(os.path.join(dirname, filename))
if os.path.splitext(filename)[1] == '.tiff':
tiff_files.append(os.path.join(dirname, filename))
for i, file in enumerate(hspy_files):
print('[', i, ']', file.split('/')[-1])
Calibrating real space pixel size¶
au_refImages = []
for file in hspy_files:
if 'au_xgrating_100kX_20cmCL_10umAp_ibf' in file:
au_refImages.append(file)
if 'au_xgrating_150kX_20cmCL_10umAp_ibf' in file:
au_refImages.append(file)
if 'au_xgrating_200kX_20cmCL_10umAp_ibf' in file:
au_refImages.append(file)
auRefIm = hs.load(au_refImages)
plt.figure()
plt.subplot(131, title = '100kX_20cmCL')
plt.imshow(auRefIm[1].data)
plt.axis('off')
plt.subplot(132, title = '150kX_20cmCL')
plt.imshow(auRefIm[2].data)
plt.axis('off')
plt.subplot(133, title = '200kX_20cmCL')
plt.imshow(auRefIm[0].data)
plt.axis('off')
The pitch of this standard cross-grating sample is known (500 nm). Here we can measure and calibrate pixel sizes - interactively:
im= auRefIm[1]
im.plot()
line = hs.roi.Line2DROI(x1=68.7, y1=213.96, x2=189.9, y2=64.167, linewidth=10)
line.add_widget(im)
trace = line.interactive(im)
trace.plot()
MAG = [
'100kX',
'150kX',
'200kX'
]
# pitch of the au cross-grating is 500 nm
nm_per_pixel = [
np.mean([500 / (101-36.5),500 / (164.7-101)]),
np.mean([500 / (109.46 - 14.91),500 / (205.97 - 109.46)]),
np.mean([500 / (134.57 - 10.31),500 / (263.27 - 134.57)])
]
FOV = np.asarray(nm_per_pixel) * 256
#dictionary of MAG cal - MAG keys and FOV in nm values
MAG_cal_FOV = dict(zip(MAG, FOV))
print(MAG_cal_FOV)
# save to a file
np.save(os.path.join(base_folder, 'MAG_cal_FOV_nm.npy'), MAG_cal_FOV)
Calibrating reciprocal space pixel size¶
Reference diffraction data, only using 20 cm CL cases:
au_refDiff = []
for file in hspy_files:
if '_sumdiff' in file:
au_refDiff.append(file)
auRefDiff = hs.load(au_refDiff)
v_max = 15
plt.figure()
plt.subplot(131, title = '100kX_20cmCL')
plt.imshow(np.log(auRefDiff[2].data), vmax = v_max)
plt.axis('off')
plt.subplot(132, title = '150kX_20cmCL')
plt.imshow(np.log(auRefDiff[1].data), vmax = v_max)
plt.axis('off')
plt.subplot(133, title = '200kX_20cmCL')
plt.imshow(np.log(auRefDiff[2].data), vmax = v_max)
plt.axis('off')
It should not matter which one we choose to analyse here. So continuing with the first pattern:
dp = auRefDiff[0]
dp = pxm.signals.ElectronDiffraction2D(dp)
dp.plot(vmax = 8000)
Lets do a rough measurement of pixel size in reciprocal space - interactively:
dp.plot(vmax=8000)
line = hs.roi.Line2DROI(x1=186.251, y1=342.793, x2=307.451, y2=193.00, linewidth=5)
line.add_widget(dp)
trace = line.interactive(dp)
trace.plot()
# centre : 240, 275.5
# 111 ring diam in pixels: 92 - radius 46
# Au 111 reflection (A^-1): 1 / 2.355
diff_cal = (1/2.355) / 46
print(np.round(diff_cal, decimals=5), 'A^-1 to pix')
The above value is a rough initial estimate. We will refine this value later in the notebook.
Diffraction roundness correction¶
First we create a smaller dataset by skipping a number of probe positions, e.g. 10 here:
file = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/expands_notebook/au_xgrating_cals/20191018 113710/au_xgrating_150kX_20cmCL_10umAp_scan_array_255by255_diff_plane_515by515_.hdf5'
dp_full = hs.load(file, lazy = True)
dp_skip10 = dp_full.inav[::10,::10]
dp_skip10.compute()
# dp_skip10.save('/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/20191018 113710/au_xgrating_150kX_20cmCL_10umAp_scan_array_255by255_diff_plane_515by515_skip10.hdf5')# dp_skip10= pxm.load_hspy('/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/au_xgrating_cals/20191018 113710/au_xgrating_150kX_20cmCL_10umAp_scan_array_255by255_diff_plane_515by515_skip10.hdf5')
Loading the mask and applying it to the skip10 data:
mask = hs.load(mask_path)
mask.plot()
dp_skip10_masked = dp_skip10 * mask
dp_skip10_masked = pxm.signals.ElectronDiffraction2D(dp_skip10_masked)
dp_skip10_masked_sum = dp_skip10_masked.sum()
We need to have the central beam at the centre of the diffraction patterns to continue with this analysis. Here we find the centre of the pattern on the sum diffraction pattern and crop around this position:
centre = pxm.utils.expt_utils.find_beam_center_interpolate(dp_skip10_masked_sum, sigma=5, upsample_factor=4, kind='linear')
print('Coordinates of the centre: ', centre)
det_size = dp_skip10_masked_sum.axes_manager[0].size
crop_size = min(int(det_size - centre[0]),int(det_size - centre[1]), centre[0], centre[1])
print('croped size of the patterns: ', int(2 * crop_size))
dp_skip10_masked_crop = dp_skip10_masked.isig[int(centre[0]-crop_size):int(centre[0]+crop_size), int(centre[1]-crop_size):int(centre[1]+crop_size)]
Now we can correct for shifts in the diffraction pattern by centring the direct beam:
dp_skip10_masked_crop.center_direct_beam(method='interpolate', sigma=5, upsample_factor=4, kind='linear')
dp_skip10_masked_crop.axes_manager[2].offset = 0
dp_skip10_masked_crop.axes_manager[3].offset = 0
dp_skip10_masked_crop.plot()
After the above correction, we can sum all the diffraction patterns and apply the reciprocal space calibration:
dp_sum = dp_skip10_masked_crop.sum()
dp_sum = pxm.signals.ElectronDiffraction2D(dp_sum)
dp_sum.set_experimental_parameters(beam_energy=300, camera_length=20.)
dp_sum.set_diffraction_calibration(calibration=diff_cal)
cal_lib = CalibrationDataLibrary(au_x_grating_dp=dp_sum)
# cal = CalibrationGenerator(calibration_data=cal_lib)
cal = CalibrationGenerator(diffraction_pattern=dp_sum)
dp_sum.plot(vmax = 300, cmap='magma_r')
Using this pattern we can correct for diffraction roundness:
cal.get_elliptical_distortion(mask_radius=10,
scale=diff_cal, amplitude=80, direct_beam_amplitude=4e3,
asymmetry=0.9,spread=2)
corr_matrix = cal.get_correction_matrix()
print(corr_matrix)
residuals = cal.get_distortion_residuals(mask_radius=10, spread=2)
residuals.plot(cmap='RdBu', vmax=0.005)
cal.plot_corrected_diffraction_pattern(vmax= 80)
Now we can refine the reciprocal space pixel size
centre = (int((dp_sum.data.shape[0]/2)), int((dp_sum.data.shape[1]/2)))
profile = epsic.radial_profile.radial_profile(dp_sum.data, centre)
fig, ax = plt.subplots(1,1)
ax.plot(profile)
ax.set_xlim([0, 50])
ax.set_ylim([0, 200])
ax.annotate('au_111', xy=(44, 95), xycoords='data',
xytext=(0.8, 0.95), textcoords='axes fraction',
arrowprops=dict(facecolor='black', shrink=0.05),
horizontalalignment='right', verticalalignment='top',
)
diff_cal_refined = (1/2.355) / 44.4
print('Initial estimate:', np.round(diff_cal,5), 'A^-1 to pix')
print('Refined value: ', np.round(diff_cal_refined,5), 'A^-1 to pix')
Scan rotation calibration¶
This calibration aims to correct for the rotation angle between the scan coil fast and slow directions and the X/Y axis of the detector. The data we load here is a 4DSTEM dataset of a largely defocused probe (i.e. Ronchigram). We then take two lines corresponding to the slow and fast scan directions.
scan_data = '/dls/science/groups/e02/Sample_data/sample_pencilbeam_SED_data/expands_notebook/au_xgrating_cals/scan_rot/20201214_153300.hdf5'
scan = hs.load(scan_data, lazy = True)
fast_dir = scan.inav[1:,60]
slow_dir = scan.inav[60,1:]
fast_dir.compute()
slow_dir.compute()
slow_dir.plot()
Cropping to the area within the bright-field disc, removing the cross, and applying a Gaussian blur to better track the motion:
index = [255,256,257,258,259]
fast_dir_no_cross = np.delete(fast_dir.data, index, axis = 1)
fast_dir_no_cross = np.delete(fast_dir_no_cross, index, axis = 2)
slow_dir_no_cross = np.delete(slow_dir.data, index, axis = 1)
slow_dir_no_cross = np.delete(slow_dir_no_cross, index, axis = 2)
fast_dir_no_cross = hs.signals.Signal2D(fast_dir_no_cross)
slow_dir_no_cross = hs.signals.Signal2D(slow_dir_no_cross)
fast_dir_crop = fast_dir_no_cross.isig[132:429, 108:415]
slow_dir_crop = slow_dir_no_cross.isig[132:429, 108:415]
from scipy.ndimage import gaussian_filter
fast_dir_crop_gb = gaussian_filter(fast_dir_crop.data, sigma = 2)
fast_dir_crop_gb = hs.signals.Signal2D(fast_dir_crop_gb)
slow_dir_crop_gb = gaussian_filter(slow_dir_crop.data, sigma = 2)
slow_dir_crop_gb = hs.signals.Signal2D(slow_dir_crop_gb)
slow_dir_crop_gb.plot()
estimate the shifts in the image:
shifts_slow = slow_dir_crop_gb.inav[:20].estimate_shift2D('current', roi = (132,292,150,300), normalize_corr=True)
shifts_fast = fast_dir_crop_gb.inav[:20].estimate_shift2D('current', roi = (132,292,150,300), normalize_corr=True)
shifts_slow = np.asarray(shifts_slow)
print(np.sum(shifts_slow, axis=0))
shifts_slow_sum = np.sum(shifts_slow, axis=0)
shifts_fast = np.asarray(shifts_fast)
print(np.sum(shifts_fast, axis=0))
shifts_fast_sum = np.sum(shifts_fast, axis=0)
plt.figure()
plt.plot([0, shifts_slow_sum[0]], [0, shifts_slow_sum[1]], '-r', label = 'slow_dir')
plt.plot([0, shifts_fast_sum[0]], [0, shifts_fast_sum[1]], '-b', label = 'fast_dir')
plt.axis('square')
plt.legend()
import scipy
rot_angle = (abs(math.degrees(math.atan(shifts_fast_sum[1]/shifts_fast_sum[0])))
+ 180 - abs(math.degrees(math.atan(shifts_fast_sum[1]/shifts_fast_sum[0])))
- abs(math.degrees(math.atan(shifts_slow_sum[1]/shifts_slow_sum[0]))))
print(rot_angle)
slow_dir_crop_gb_rot = scipy.ndimage.rotate(slow_dir_crop_gb.data, angle=rot_angle, axes = (1,2))
slow_dir_crop_gb_rot = hs.signals.Signal2D(slow_dir_crop_gb_rot)
slow_dir_crop_gb_rot.plot()
fast_dir_crop_gb_rot = scipy.ndimage.rotate(fast_dir_crop_gb.data, angle=rot_angle, axes = (1,2))
fast_dir_crop_gb_rot = hs.signals.Signal2D(fast_dir_crop_gb_rot)
fast_dir_crop_gb_rot.plot()
With the parameters reached in this notebook we can now calibrate all the subsequent experimental SEND datasets.